Density Fluctuations in an Electrolyte from Generalized Debye-Hiickel Theory 
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Near-critical thermodynamics in the hard-sphere (1,1) electrolyte is well described, at a classical 
level, by Debye-Hiickel (DH) theory with (+, — ) ion pairing and dipolar-pair- ionic- fluid coupling. 
But DH-based theories do not address density fluctuations. Here density correlations are obtained 
by functional differentiation of DH theory generalized to non-uniform densities of various species. 
The correlation length £ diverges universally at low density p as (Tp) -1 ^ 4 (correcting GMSA theory). 
When p = p c one has £ ~ CtTA 1 ^ 2 &s t = (T — T c )/T c — » 0+ where the amplitudes fo~ compare 
informatively with experimental data. 
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There is a major puzzle in the theory of fluid crit- 
icality in model ionic systems jl| because experiments 
jl],D reveal that certain electrolyte solutions exhibit clas- 
sical (i.e., van der Waals as against the usual Ising-type) 
critical exponents over as much as 1 to 3 decades when 
\t\ = \T — T c \/T c — > 0. Probably there is always a scale 
t x below which the behavior crosses over from classical to 
Ising; but attempts to explain how t x might vary from 
~ 1 to ~ 10 have so far been unconvincing. Initial 
efforts have addressed the simplest case: the 'Restricted 
Primitive Model,' consisting of N = N + + N~ = V p 
hard-sphere ions of diameter a, carrying charges ±q in 
a medium of dielectric constant D. The hope has been 
to decide the universality class (and crossover scale t x if 
appropriate) of the RPM |g. 

To that end Fisher and Levin j|] have shown that the 
original Debye-Hiickel (DH) theory Q provides a remark- 
ably good, albeit classical account of the critical thermo- 
dynamics as judged by current simulations How- 
ever, for a satisfactory description, pure DH theory must 
be extended (i), following Bjerrum [Q,[j|, by recognizing 
bound, neutral but dipolar (+,-) ion pairs in equilibrium 
with the free ions, (ii) by including the dipolar-ionic (DI) 
solvation free energy $W , and (iii) by allowing for hard- 
core (HC) repulsions. In terms of T* = ksTDa) 'q 2 , these 
DHBjDIHC theories yield critical points in the range 
T* ~ 0.052 to 0.057 as compared with 0.052-0.056 from 
recent simulations [3(b)]. 

Now, following Ebeling and Grigo Q , one can also pur- 
sue theories based on the mean spherical approximation 
(MSA); but, for reasons currently obscure, such theories, 
even when improved in various ways [3(b), 5-8] yield esti- 
mates for T* too high by 35-50% H. Note also that the 
hypernetted chain (HNC) and other integral equations 
have no solutions in the critical region! Further study of 
the DH-based theories is thus well justified. 

To understand properly the nature of a critical point 
one must go beyond thermodynamics to study the order- 
parameter fluctuations. But, even for ionic fluids, the 
order parameter is just the overall density. Now DH- 
based theories illuminate the Debye-screening of the bare 
Coulombic potential, as seen in the exponential decay of 



the charge-charge correlations, but, unfortunately, they 
say little about the overall density- density correlation 
function, G pp {r) = (p(r)p(O)) - p 2 = p{S{r) + p[g 2 (r) - 
1]}. Our aim here is to rectify this deficiency. 

Note, especially, that the Fourier transform of G pp (r) 
yields the k-dependent susceptibility 



X (k) = G w (k)/p = X (0)/[l + ^ + 



(1) 



which diverges at k = at criticality. Indeed, xQk) de- 
termines the critical opalescence and turbidity [0 and 
specifies the (second-moment) correlation length £(T, p) 
which diverges as £g~ ft" when t — ► 0+ at p = p c ; fur- 
thermore, x(k) approaches the reduced compressibility 
y(0) = pksTKT (or its solution analog) when k — > 

As stressed by Fisher and Levin [3 (a), 7], it is valu- 
able to know the critical amplitude £q" even within a 
classical approximation since, via the Ginzburg criterion, 
that offers a possible route for estimating a crossover 
range, ±t x , outside which closc-to-classical critical be- 
havior might be realized |TT[| . 

In this Letter we show how DH theory can be gener- 
alized so as to yield, in a natural way, density fluctu- 
ations diverging at criticality Jl2|-[l4|]. The method ex- 
tends straightforwardly to the full DHBjDIHC theories 
H as we illustrate below. In particular, we calculate the 
correlation length £(T, p) explicitly within the simplest 
generalized (GDH) theory, and numerically, at improved 
levels of approximation. At low densities a novel, uni- 
versal divergence of £(T,p) is predicted for all T. In the 
critical region the results arc, as expected, classical with 
v = h; but the amplitudes £q prove informative and 
are compared with experimental data [p| in classical and 
Ising domains pd| |. 

Explicitly we proceed, following ||, by approximat- 
ing the total Helmholtz free energy F(T, p) by a sum of 
terms representing ideal gas, ionic fluid, dipole-ion, and 
hard-core contributions; but we now aim for a functional 
j3F[{pj}] = J d d rT where pi(r) = p+(r) + P-(r) and 
P2(r) are slowly varying local free-ion and dipolar pair 
densities, while j3 — 1/kgT. Since we wish to probe 
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only the density fluctuations, we follow DH theory and 
maintain electroneutrality, p+ (r) = p_ (r) , on long length 
scales. Of central concern are the (reduced) direct corre- 
lation functions which follow by functional differentiation 
(with i,j — 1,2) as 



Cij(r- 



6 2 (3F/5 Pi (r)6p 3 (r')\ 



PaO")=P*' 



(2) 



where the p\ (A = +,—,2) are the overall equilibrium 
densities. Note that the various terms in T contribute 
linearly to the Cy and, in particular, Cy (k) = <5y / 'pj. 
However, since the total local density is p(r) — pi(r) + 
2p 2 (r) |lj| one finds, with the aid of the Ornstein-Zernike 
(OZ) matrix relation for Cy(k) pC(| , the result 



Cn(k)C 22 (k)- [Ci 2 (k)] ; 



PX(k) G pp {k) 4C u (k) - 4C 12 (k) + C 22 (k) : 



(3) 



from which £ follows by expansion in k. More expedi- 
tiously one may impose infinitesimal density variations 



Pj (R) = pj [1 + Aj cos k • r] , 



(4) 



and expand the reduced free-energy density (3F/V 
in powers of the Aj: the quadratic term is then 
\ Si j PiPjCij{k)AiAj, from which the Cy follow. 

Evidently, the crucial issue is to extend DH theory to 
nonuniform but slowly varying mean densities of the var- 
ious species. Note first that the free-ion contribution be- 
comes Jlq] , via the Debye charging process fl , 



F 



■DH 



cfVipi(ri) / d<7iV>i(ri,9i), 



(5) 



where (ri; gi) is the mean electrostatic potential at the 
site ri of a fixed ion due to all the other ions when each 
carries charges ±q. If 0(r,ri) is the mean electrostatic 
potential at a general point r when the ion 1 is fixed at 
r i; one has [§ ^i(ri) = lim r ^ ri [0(r; r x ) - qx/D\r - r x |]. 
DH derived their celebrated equation for </> by approxi- 
mating the probability density for a particle of species 
A (= +,-,2) with charge q x by p x exp[- pq\(j>(r)] ~ 
Pa[1 — /3(7A0(r)] |§. (Note g 2 = 0.) In the same spirit we 
now propose to replace the constant partial (species) den- 
sity p\ by p\ (r) the (slowly varying) non-uniform density 
p6| . Our generalized (GDH) equation then reads 

[V 2 - « 2 (r)0i(r)]^(r; r x ) - -47rg x *(r - n)/D, (6) 

where, utilizing d(y), the Heaviside step function, 
0i (r) = 9(\r — ri\ — a) embodies the crucial hard-core 
boundary condition [Q, while the spatially varying coef- 
ficient 

k 2 [{p 3 }] = 47r/?£ A ai PA (r)AD = 4n(3q 2 Pl (r)/D, (7) 

reduces to the standard expression for k 2 , the inverse 
Debye length squared, when p\(r) = p\ is constant B. 



To solve (q), we adopt (Q) and expand <j> in powers of 
Ai. The coefficient </> n (r, ri) of A™ can be computed re- 
cursively, setting k = k and ri = 0, with the aid of the 
Green's function 



OO , i \ 



(8) 



where, employing modified spherical Bessel functions, 



&(a>s')=-w ~ 



s, s < x, 



=(2£+l)- 



where x = «a, s> 
Gt(s,a') 



2£+l 



s< < x < s>, (9) 
= min(s, s'), while, 
fei(*)^(« / ) + i/(s<)^(« > ), (10) 



max(s, s'), s< 
i«fi fa) 



for s, s' > x, and Pi{p) denotes a Legendre polynomial. 

Substituting in (||) and expanding C P(9 = 1/G pp to 
0(k 2 ) yields £ 2 . This requires only the £ = and 1 
terms in (||). Consequently, within pure DH theory, the 
simplest level of approximation, we find the correlation 
length is given explicitly by (recalling x = na) 



24T*x 2 



In 



(l + x) 
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x — 5x 



(l + x + ±x 2 ) 9 



2(1 



(11) 



where l/x DH (0) = 1 -x/4T*(l + x) 2 @. 

Now corrections to this result enter only in 0(x 2 ) = 
O(p), i.e., beyond the leading low-density behavior 
which, in fact, exhibits the novel divergence 



f (T, p) = l(b/36np)^{l + \nb + 0(p*)], 



(12) 



when p — » 0, where 6 = q 2 / DIzbT is Bjerrum's length. 
This expression for the density-density correlation length 
is independent of the hard-core diameter a and is thus 
universal^. We believe it represents the exact limiting 
behavior not previously noted. At low densities Debye's 
screening length £d = 1/k controls the decay of the 
charge correlations [4(b), 17]. It also diverges universally 
when p — > 0; but since we find £ w x/b^o/48, the den- 
sity correlations then decay on a shorter scale than the 
charge correlations. 

Our conclusion ( p2] ) can be checked further by using 
the HNC relation cy m -/?«y + ^h 2 j [4(b), 17,18], which 
is probably generally valid in the low density limit p^ ] 
when hij = gij — 1 decays e 
here. This leads to ||Q 

l/x(k) w 1 - i^btan-^fc^^/fc (p->0), (13) 

which, on expansion to order k 2 , yields the DH limit- 
ing law for x(0) and reproduces (p2|)! However, the true 
correlation length, ^(T, p), that determines the OZ-likc 



when h^ = g^ — 1 decays exponentially fast, as expected 
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exponential d ecay of G pp (r) is determined by the domi- 
nant zeros of (|13]). These give the different expression 



^ D {l + 2exp[-4/(^ 3 ) 1 /2] + ...} ; 



(14) 



than e — r /r , the squared charge-charge correlation. 

By contrast to (|l2| ) and (|l3j), one finds Jl6| that 
the GMSA or generalized mean-spherical approximation 
@,H predicts l/%(k) « 1 - ±Ko/[2 + fc 2 a/K] when p -> 0. 
This gives x(0) correctly to (^(p 1 / 2 ) but leads to 

£gmsa « CooGMSA « (±o&>) 1/2 = \{c?/*pb)V\ (15) 

Thus £gmsa also diverges as p -1 / 4 , but the power of 
T differs and the result is non-universal, depending on 
a. This reveals an unsuspected defect of the GMSA 
p0| . which was especially devised to satisfy a variety of 
correlation-function sum rules fL8fl . (The original MSA 
gives only a hard-sphere result for G pp (r): see, e.g., [^0[.) 

In the critical region the pure GDH result ( p"i"| ) diverges 
with exponent 2v = 1 at T* = jg, x c = 1, p* = 1/QAtt. 
The correlation length amplitude is found to be 

(£o + /a) dh = [1 + f In 2 - 6 In |] 1/2 ~ 0.7329. (16) 

This is surprising ly close to the GMSA value 0.75 [|, 
although T* and p* differ significantly § . 

However, although the pure GDH theory based on (||)- 
(0) is sufficient at low density, one must, as mentioned 
in the introduction include ion pairing to study the 
critical region. Bjerrum's ansatz for the association con- 
stant is appealing but Ebeling's result is superior j^] and 
used here. (Near criticality the numerical changes are 
minor.) In simple "DHBj" theory the ion pairs are sup- 
posed ideal |3| and one finds C22 = I//O2, C12 = 0, and 
(7xi (k) is unaltered. But that is too naive and proves un- 
physical: it is essential to include the dipole-ionic fluid 
(DI) interactions |J. 

We calculate the new non-uniform DI effects by us- 
ing the GDH equation, (^]), but with a dipolar source 
term, i.e., + and — point charges at ri = ±iai, where 
ai specifies the orientation and typical charge separa- 
tion, ai(T) = |ai| ||. The associated bispherical ex- 
clusion zone is approximated by a sphere of radius ai 
||. Thus the Green's function (||) still applies, but with 
x — > x-2 = Ka2- At low T, ai = a ("contact'') and 
02 = 1.16198a (angular average) are reasonable || and 
the sensitivity to these values is readily tested. 

Angular integration over the dipole orientations 
is complicated, yielding the solution 4>dip( r ', r i) a s 
a multiple sum with Clebsch-Gordan coefficients, 
Ce 1 .i 2 (mi , rri2 \l , m). To obtain V>2 (r; q) for use in the pair 
analog of (j^), the self-potential of the source dipole is 
subtracted. To order k 2 one needs only I = 0, 1, 2, which 
give explicit expressions jL6[ with low-density expansions 



Cf/(k) 



1 + 



[16,19] which diverges as (T/p) 1 / 2 . Thus G pp (r) has a 
small but longer-range tail decaying slightly more slowly ~ DI 



C12 (k) 



Xx\p 2 

20T*p 2 
k 2 

36 K 2 



k 2 



Mxo 
27 J-2 



25^,2 
21 X 2 



15 X 2 



Q(x 



-0{xi) 
0(fc 4 )}, 



(17) 



x 2 - |a?| 



|x 2 



1-3 
5 X 2 



■§ x 2 



Q(x 



0(* 4 )] + 0(fc 4 )}, (18) 



where x\ = na\ while C^ 1 (k) = 0. 

Finally, hard-core exclusion may be approximated by 
local, free- volume terms T HC = — X)j Pi m [l — Tlj BjPj] 
with (i) B\ = \B2 — 4a 3 /3v3 to yield bcc close-packing 
or (ii) B\ = \f$2 = 27ra 3 /3 for the exact ion- ion second 
virial coefficient [||; or (iii), perhaps preferably, by the 
local Carnahan-Starling mixture form plfl 



(Co - J) ln(l - Cs) 



CI 



3CiC2 
l-Ca Ca(l - Ca) 2 



(19) 



where ( n = Pi{ r ) a i with o"j the hard-core diam- 

eter of species i; we take af — \a\ 
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For densities 



near critical, only the second virial coefficients prove sig- 
nificant. Being local, the approximations (i)-(iii) give 
Cfj C '(k) independent of k. Nonlocal effects are easily 
included at the second-virial-coefficient level; but that 
changes the critical amplitude £q~ by less than 1%. 
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FIG. 1. Inverse square of the density-density correla- 
tion length, £(T, p), on the critical isochore according to: 
(a) pure DH theory; (b) with Bjerrum association and 
dipole-pair-ionic-fluid coupling (DHBjDI) with ai = a, 
a2 = 1.16198a; (c)-(e) with hard-core terms (i)-(iii): see text. 

For the DHBjDI theories the equilibrium equations re- 
quire numerical solution. Fig. 1 shows the resulting in- 
verse square correlation lengths vs. T* on the critical 
isochore for various levels of approximation. The linear 
approach of all plots to £~ 2 = represents the expected 
classical prediction v = ^ . The effects of the DI coupling 
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are less dramatic near T c than might have been guessed. 
With the assignments a\ = a, a 2 — 1.162a || the crit- 
ical amplitudes are £/a ~ 0.7511, 0.7776, 0.8186, and 
0.8147 for pure DHBjDI theory and with HC treatments 
(i)-(iii), respectively. Increasing a\ to 1.15a lowers £g by 
no more than 8.3%. Similarly, taking a-i to be 1.150a 
leads to a reduction of less than 1.1%. (The correspond- 
ing changes in T* , p*, etc. can be found in [3(b)].) 

To compare our results for £q" with experiments on 
systems that might plausibly be modeled by the RPM, 
one needs not only data for Q but also some esti- 
mate of the effective hard-core diameter, a. That might 
be obtained by matching the p* predictions to experi- 
ment. To that end, we re-express our results above as 
£+Pe /3 ~ 0.2275, 0.2302, 0.2375, and 0.2368 (in contrast 
to 0.1251 for pure DH theory and 0.1828 for the GMSA 

0)- 

Beyond that one must recall that our theory is clas- 
sical with v = vm f = 5 whereas observations indi- 
cate the Ising value vj s ~ 0.63 or crossover to that for 
t < t x . However, for 3-D lattice gases, which are de- 
scribed by vis for all t < 1, the mean-field estimates 
for £q , say £q IF , agree with numerical estimates, say 
£q s , to within 10% |2^] (for various lattice structures). 
On the other hand, if crossover is found, the fits for 
t < t x and t > t x should roughly obey the match- 
ing formula ^ S /^ F = {tx) Vls ~ UMF ■ Data for Na+ND 3 
with t x — (7-9) x 10 -3 have been fitted in both regions 
[2(a) (c)] and confirm the relation. For this system we 

find £o s Pc /3 = 0.34 ± 0.03, which is some 40 to 50% 
above our estimates. For Pitzer's salt [2(b)(e)] we may 
postulate a crossover at t x — 1 x 10~ 4 M: this yields 

£o s Pc = 0.21 ± 0.04 which encompasses our values. 
Tetra-n-butyl-ammonium picratc in n-tridecanol [2(d)] 

displays crossover and we find £,Q S pl^ 3 — 0.22, close to 
our prediction. Finally, for the same salt in other sol- 
vents [2(f)] quite similar values of £q s fit the turbidity 
data. Overall the agreement is encouraging when using 
the Ising-fitted amplitudes. It is puzzling, but perhaps 
significant, that in the mean- field region outside t x the 
values of £q IF pV Z are all larger (by factors of 2-3) than 
our classical theory predicts! 

In conclusion, we have shown how to calculate density- 
density correlations within DH theory and its extensions 
j|. At low densities the correlation length, £(T,p), di- 
verges in unexpected but universal fashion potentially 
amenable to experimental check. In the critical region 
comparison with experiments on electrolytes proves in- 
structive and raises further questions. More concretely, 
the present theory enables the observed classical-to-Ising 
crossover to be addressed via the Ginzburg criterion 

10- 

A naive extension of the functional approach out- 
lined is not sufficient for studying the charge correlations 
at higher densities since, when k — * 0, the associated 
density perturbations, even when infinitesimal, induce 



a macroscopic charge imbalance. However, approaches 
which separate out the long-range Coulombic contribu- 
tions should lead to progress. 
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